1/f noise and multifractality in atmospheric-C02 records 
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We study the fluctuations in the measured atmospheric CO2 records from several stations and 
show that it displays 1/f noise and multifractality. Using detrended fluctuation analysis and wavelet 
based methods, we estimate the scaling exponents at various time scales. We also simulate CO2 
time series from an atmospheric chemistry-transport model (CTM) and show that eventhough the 
model results are in broad agreement with the measured exponents there are still some discrepancies 
between them. The implications for sources and sinks inversion of atmospheric-C02 is discussed. 
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Atmospheric carbon dioxide is thought to be one of 
the greenhouse gases primarily responsible for the current 
phase of global warming 0. In view of this significance 
for the global climate and life on earth, considerable ef- 
fort has gone into measuring CO2 from various platforms 
in the last several decades. Direct record of atmospheric 
CO2 is available since the 1950s and hourly time series are 
constructed more recently at several measurement sta- 
tions In the context of atmospheric systems, where 
the mean behaviour is reasonably well known from such 
long term records, the fluctuations decide the actual out- 
come. However, the fluctuation properties of CO2 time 
series have not been studied in sufficient detail. Such a 
study would help identify the class of process, e.g, corre- 
lated random walk, Levy flight etc., underlying the CO2 
biogeochemistry and will provide a clue if it is a self- 
organised critical system. Secondly, this would serve as 
a diagnostic tool to compare the outputs from ab initio 
atmospheric CTMs as well as to constrain the processes 
that could be represented in them. This could lead to 
better understanding of the global CO2 cycle. 

Scale invariance and 1// noise are considered to be 
the signatures of complex systems In recent years, 
several time series from physiological, economic, physi- 
cal systems including those of atmospheric parameters 
have been subjected to scaling and spectral analysis 
If x(t) be the given time series, then the scaling can be 
expressed as, x(at) k a H x(t), where the symbol w de- 
notes statistical equality and the the scaling exponent H 
is also called the Hurst exponent. In many of naturally 
occuring phenomena, a suitable measure of their fluctu- 
ations F(L) displays a power law as a function of length 
L of the time series, i.e, F(L) oc L H . Physically this 
implies that if < H < 1 then they are long range corre- 
lated. In particular, H =1/2 corresponds to white noise 
and H = 1 to 1/f noise. A suitable stochastic model 
of long range correlations is the fractional Brownian mo- 
tion (fBm) 5] . Significantly, many natural processes such 
as the temperature and precipitation (as also many eco- 
nomic and physiological time series) behave like a fBm 



process and exhibit long range correlation with a expo- 
nent that is claimed to be universal in some of those cases. 
For instance, in the case of temperature fluctuations at 
inland sites far from the oceans the exponent is H « 0.65 
0. See also Ref Q for a debate. Similar analysis of 
rainfall has been used to suggest that rain events might 
be a case of self-organised criticality (SOC) Infact, 
the idea that many natural processes might display 1/f 
noise and SOC-like behaviour is one of the impetus for 
this line of research. However, to the best of our knowl- 
edge, the fluctuation properties of the chemical species 
in the atmosphere such as CO2 have not been studied. 
In this letter, we examine the experimentally measured 
high frequency (hourly) CO2 time series from 28 stations 
and show that their flucuations display a power law with 
H « 1 corresponding to / _/3 (j3 w 1) noise at interme- 
diate frequencies. We also simultaneously examine the 
model simulations of atmospheric-C02 and discuss its 
implications. 

The measured CO2 concentrations constitute a non- 
stationary time series and hence we apply the detrended 
fluctuation analysis (DFA) [lj} to quantify the fluctua- 
tions and its behaviour in the spectral domain. An al- 
ternative approach is to use wavelets [Tlj of appropriate 
basis such that it allows for separation between high and 
low frequency components through successive levels of 
decomposition. We will only briefly mention the tech- 
niques here since both the methods are widely reported 
in the literature 0, 0. We denote by x{k) the CO2 
concentration in parts per million (ppm) measured at 
discrete times k = 1,2, ....k max . The corresponding inte- 
grated series given by, y(t) = Y%=1 t = !> 2 -- fc m«i- 
The detrending is done by dividing the series y(t) into M 
segments each containing r data points and by empiri- 
cally fitting a linear, quadratic, cubic ... function y(t) 
to each of them. Then, the fluctuations about the fitted 
trend in each of the segment is obtained as, 
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FIG. 1: (a) Hourly time series of measured CO2 concentra- 
tion at Brotjacklriegel, Germany, (b) DFA applied to this 
time series (solid dots) shows two distinct time scales with 
different scaling exponents. The blue line is the best fit line 
with slope of 1.440 ± 0.012. The black line fits the second 
region with slope of 0.953 ± 0.011. The dashed line is the 
shuffled data and has a slope of 0.517, close to that of white 
noise, (c) The power spectrum of CO2 series displays corre- 
sponding distinct frequency scaling. The blue line is obtained 
as the best fit one with slope -1.920 ± 0.023 and the black 
line has the slope of -1.080 ± 0.032. The latter region dis- 
plays low-frequency 1/f noise, (d) Plot of H q against q as a 
signature of multifractality. 



The fluctuations F(t) averaged over all the M seg- 
ments of size r is denoted by (F(t)). In many cases, 
(F(t)) oc t h . In the case of standard random walk, 
H = 3/2. Note that the exponent H is related to the 
power spectral exponent of x(k) through (3 = 2H — 1 and 
to its fractal dimension D = 2 — H [jj. We also obtain the 
fluctuation function F(t) using Daubechies wavelet [ll| . 
an orthogonal wavelet family with largest number of van- 
ishing wavelet coefficients. We perform 9 levels of wavelet 
decomposition corresponding to quadratic detrending. 

The hourly CO2 data are available in the public do- 
main for 28 stations, some of which have long period of 
measurements going as far back as 1972 |3j. As a typical 
case, Fig ^a) displays the time series of CO2 concen- 
trations in ppm obtained at Brotjacklriegel, Germany. 
Apart from the hourly fluctuations, the profile shows an 
approximate periodicity of about one year which is re- 
lated to the consumption pattern of CO2 by the vegeta- 
tion in the northern hemisphere. In Fig^b), we show 
the result of DFA performed on this data. Two distinct 
scaling regimes can be identified. For a time scale of few 
hours to two weeks, denoted by blue line in Fig^b), 



we obtain a best fit slope of 1.440 and this closely corre- 
sponds to 3/2 one gets for a random walk process. For a 
time scale of about 30 days (shown as black line), the esti- 
mated slope is 0.953 and this regime closely corresponds 
to 1/f noise. At longer time scales beyond t > 300 
days, we see an approximate white noise like behaviour 
with H « 0.512. In order verify that these correlations 
are genuine, we create a surrogate by randomly shuffling 
the series while retaining the same distribution. Ran- 
dom shuffling would kill the correlations present in the 
data. The DFA performed on this shuffled surrogate data 
(shown as a dashed line) has a slope of 0.517, in close 
agreement with the white noise exponent of 0.5. We 
have also verified the DFA exponents using the wavelet 
decomposition as well. The power spectrum in Fig^c) 
displays, as expected based on DFA results, 1/f noise at 
intermediate frequencies (black line) and 1 / f 2 behaviour 
at higher frequencies (blue line). Notice that the esti- 
mated exponents H in Fig^b) are related to the power 
spectral exponents (3 in Fig IHc) through j3 = 2H — 1 . 
The results in Fig^c) also indicate that at low frequen- 
cies the spectrum is approximately flat representative of 
white noise. This succession of three different regimes 
in the power spectrum corresponding to white noise (/°) 
at low, 1/f noise (/ _1 ) at intermediate and random walk 
(J -2 ) at high frequencies can be generated from a super- 
position of exponential relaxation processes of the type 
N(t) — iVoexp(— Xt) with the uniform distribution g(X) 
of relaxation rates A 01 ■ While it might be tempting 
to suggest this as a model for atmospheric-C02 variabil- 
ity, further work needs to be done both at the level of 
statistics and much more in terms of CO2 science to un- 
derstand the modeling aspect further. Further, we will 
present results to show that this scaling behaviour is 
generic to hourly CO2 data from all the stations, with 
the notable exception of Samoa and South Pole. 

As shown above, Fig^b) displays atleast three distinct 
scaling regimes which indicates the possibility that CO2 
time series could be a multifractal. In a monofractal, all 
the moments of the fluctuations have the same scaling 
exponent H. For a multifractal, qth moment is given by, 
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and the exponents H q are dependent on order q. We com- 
pute H q using multifractal-DFA (MFDFA) and wavelet 



formalism [13j and in Fig IHd), we show H q plotted 
against q. Note that at q = 2, we have H2 — 1.062, in 
close agreement with H obtained using DFA and wavelet 
methods earlier. We will show that CO2 series from all 
the stations display multifractal characteristics. 

Our aim in rest of the study is to show that the scal- 
ing results presented in Fig. ^ are generic to almost 
all the sites for which the data is available. For this 



purpose, we consider the following 8 stations [T^ cho- 
sen to represent various geophysical features; namely, 
Waldhof (LGB), Schauinsland (SSL) both in Germany, 
Ryori (RYO), Hateruma (HAT) both in Japan, Barrow 
(BRW) in the USA, Jubany ( JBN) in Argentina, Puszeza 
Borecka (DIG) in Poland and Anmyeon-do (AMY) in Re- 
public of Korea. We show the DFA fluctuation functions 
for these stations in Fig|2ta). In all the cases, atleast 
two distinct scaling regimes are visible; the blue line 
represents the random walk type scaling regime while 
the black line represents the 1// noise type scaling. 
The mean dfa exponent for these stations in the former 
regime is H w 1.404 ± 0.125 while for the latter it is 
H w 0.978 ± 0.087 closely approximating unity neces- 
sary for l/f noise. The power spectrum shown in Fig 
Efb) provides a direct evidence of this. Beyond this, at 
longer time scales of about 130 days (logr > 3.5), some 
of the stations, AMY for instance, display white noise 
type scaling as indicated by the dashed line. However, 
the white noise regime does not seem to exist in all the 
stations. In order to obtain a global picture, the his- 
togram of DFA exponents for all the stations is shown in 
Fig [3] Note that most of the exponents assume a value 
between 0.8 and 1.1 which is our central result that in- 
termediate frequency CO2 series exhibits l/f noise. The 
power spectrum of these stations (not shown here) also 
confirms the l/f scaling. In Fig[3Jb), the histogram of H 
for short time scales of less than 14 days (l/f 2 regime) 
is shown. Most of the exponents lie within the range 1.3- 
1.7, in close vicinity of 1.5 needed for l/f type power 
spectrum. 

Thus, Figs and taken together reveal l/f noise 
(indicated as black line) for a time scale of approximately 
30 - 50 days in high frequency CO2 records and it joins 
the host of other phenomena that display l/f behaviour 
0. Performing similar analysis on the daily averaged 
data (not shown here) shows up l/f noise on similar 
time scales. 

The theoretical modelling of CO2 fluxes is presently 
an important activity and it reflects the current under- 
standing of global CO2 cycle 0|. In order to compare 
our observation based results against a theoretical bench- 
mark, we perform scaling analysis using DFA on the 
time series obtained from atmospheric CTM. We sim- 
ulate the atmospheric-C02 time series using a transport 
model and prescribed surface fluxes due to oceanic ex- 
change (monthly), terrestrial biosphere (3-hourly) and 
annual mean fossil fuel emission [lfij ]. The transport 
model is based on the CCSR/NIES/FRCGC atmospheric 
general circulation model (AGCM) nudgged with NCEP- 
FNL analysed wind vectors and temperature [l7|. The 
AGCM includes transport of chemical constituents due 
to advection, convection, and turbulent mixing processes. 
The model is run at T106 (~1. 125x1.125) horizontal res- 
olutions and 32 vertical layers for the period 2002-2003 
following a spin-up of 2 years. The DFA exponents ob- 




FIG. 2: (a) Fluctuation function F(r) shown in log-log plot 
for seven stations. Atleast two distinct scaling regions can be 
seen. Blue line has a slope of 1.5 and the black line has a slope 
of 1 meant for comparison purposes, (b) Power spectrum for 
the same set of stations as in (a) in the same order in which 
the station name appears in the legend, y-axis values are 
shifted for clarity. The slope of black line is -1 and that of 
blue line is -2, both meant for comparison with the power 
spectrum of actual data. 
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FIG. 3: Histogram of DFA exponents for all the 28 stations, 
using both model simulations (blue) and observations (ma- 
genta), (a) shows exponents from l/f region and (b) shows 
exponents from l/f 2 exponents (see text for details). 



tained from the CO2 time series simulated by the model 
at hourly time interval are shown in Fig. [3] Most of 
the model simulated exponents lie within the bins 0.7-1.2 
(shown as blue line) though about half of them assume 
values in the range of 0.7 to 0.8. The overlap between 
the two histograms in Fig a) shows an overall broad 
agreement along with the differences in details that re- 
flects partly the modelling deficiencies. For instance, it is 
known that the correlation between the measured CO2 
fluxes and the simulated ones for seasonal time scales 
(corresponding partly to l/f scaling regime) is 0.8 or 
greater |l^. The histograms for the f~ 2 random walk 
type regime in Fig |3{b) shows that even more discrep- 
ancies exist between measured and model results. This 
is because the synoptic scale (~7-10 days) variations in 
CO2 are not reproduced by the model simulations (cor- 
relation 0.5 or smaller) 15]. With this short period of 
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FIG. 4: Dimension spectrum as a signature of multifractality 
at intermediate frequencies (1// regime at q — 2) for the same 
set of stations shown in Fig [3 At g = 2, most of the stations 
have H2 » 1. 



simulation experiment, namely 2 years, it is not possible 
to observe the white noise regime. 

The results presented above are generic for most of 
the stations except remote sites like Samoa and South 
Pole. Significant deviations are observed mainly in the 
f~ 2 regime at remote observation sites, roughly in pro- 
portion to the distance from the regions of fastest CO2 
exchange, namely the forested lands and countries with 
greater human activity. For instance, the large devia- 
tions between the scaling exponents from modelled and 
observed time series are seen for the periods until 2, 7, 
8 and 74 days at Izana, Samoa, Mauna Loa and South 
Pole, respectively. 

This comparison results suggest that the transport and 
flux models are fairly realistic over the region of fast and 
large CO2 exchanges at the surface. But further studies 
are required to understand the behaviour of CO2 varia- 
tions at the remote sites like Samoa and South Pole. Both 
these stations display significant deviations from the scal- 
ing results presented above. We believe, such information 
can be used as a diagnostic tool for testing the terrestrial 
and ocean ecosystem flux model results as appropriate. 
In sources and sinks inversion of atmospheric-C02, the 
forward model simulation lengths of basis functions and 
presubtracted fluxes are dictated by the flux memory in 
CO2 time series Thus such scaling analysis would 

help in defining the forward simulation lengths and the 
time scale of future flux inversion studies. 

Further, we show that multifractality is generic to CO2 
time series. We calculate H q for the observed data at 
stations shown in Fig fusing MFDFA and wavelet for- 
malism. We stress that the results from multifractal-DFA 
[l^j l and wavelets are quantitatively similar. The MFDFA 
results from modelled CO2 time series also exhibit mul- 
tifractality, though they do not agree quantitatively with 
the observation based results. This once again suggests 
that the models are picking up broadly the correct be- 
haviour in atmospheric-C02 but missing out the details. 

We have studied the scaling properties of high- 
frequency atmospheric-C02 time series at 28 stations 



around the globe. At intermediate time scales of about a 
month, CO2 data exhibits 1/f noise. At shorter time 
scales of about a week, we get 1/f 2 noise. We also 
show that the CO2 series is a multifractal implying that 
manty different time scales are present in the system. We 
also examine the ab initio atmospheric chemistry model 
outputs and show that the model results reasonably re- 
produce the statistical properties of the measured series. 
The implications of these results in estimation and anal- 
ysis of surface fluxes are discussed. 

This work is possible due to high precision measure- 
ments of atmospheric-C02 by several groups around the 
globe (see the full list in WMO/GAW/WDCGG, 2006). 
Their relentless effort is greatly appreciated. We thank 
Hajime Akimoto for supporting this research. 
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